December 16, 2022

Outline

  • Features

  • Overview

  • Toy example

  • Actual case study: The Mitchell River, Australia

  • Shiny Application

Why ?

Features


  • Powered by Rcpp and RcppArmadillo

  • Find optimal solutions

  • Solve a specific new problem using mathematical programming

  • Ideally, fast 🙂

Overview

Three questions: Planning objective, Optimization objective, and Intensity of threats.

Workflow

Q1. Planning objective?

Recovery vs Conservation

Workflow

Q2. Threat intensities?

Binary vs Continuous


Bolch, E.A. et al. (2020). Remote Detection of Invasive Alien Species

Q2. Threat intensities?

Workflow

I) Input data validation

Inputs (inputData())

Data that defines the spatial prioritization problem (planning units data, feature data, threats data, and their spatial distributions).

5 + 2 (optional)
inputData(
  pu,
  features,
  dist_features,
  threats,
  dist_threats,
  sensitivity = NULL,
  boundary = NULL
)

Inputs

What if I have my data in another format?
A: There are functions that allow you to import and convert data to data.frame format

Examples:

##from .txt to data.frame
df <- read.table('myfile.txt',sep='\t')
##from .csv to data.frame
df <- read.csv("myfile.csv")

Source: https://libguides.umn.edu


#1 Planning units data (pu)

data.frame


Specifies the planning units (PU) of the corresponding instance and their corresponding monitoring cost and status. Each row corresponds to a different planning unit. This file is inherited from the pu.dat in Marxan.

  • id
    integer unique identifier for each planning unit.

  • monitoring_cost
    numeric cost of including each planning unit in the reserve system.

  • status (optional)
    integer value that indicate if each planning unit should be available to be selected (0), locked-in (2) as part of the solution, or locked-out (3) and excluded from the solution.

## id monitoring_cost status
## 1                2      0
## 2                2      0
## 3                2      0
## 4                2      0
## 5                2      0
## 6                2      0

#1 Planning units data (pu)

id


Game, E. T. y H. S. Grantham. (2008). Manual del Usuario de Marxan: Para la versión Marxan 1.8.10

#1 Planning units data (pu)

#1 Planning units data (pu)

#2 Features data (features)

data.frame


Specifies the conservation features to consider in the optimization problem. Each row corresponds to a different feature. This file is inherited from the spec.dat in Marxan.

  • id
    integer unique identifier for each conservation feature.

  • target_recovery
    numeric amount of recovery target to achieve for each conservation feature.

  • target_conservation (optional)
    numeric amount of conservation target to achieve for each conservation feature.

  • name (optional)
    character name for each conservation feature.

##id target_recovery     name
## 1              11 feature1
## 2              16 feature2
## 3               8 feature3
## 4               9 feature4

#2 Features data (features)

#2 Features data (features)

#2 Features data (features)

#2 Features data (features)

recovery_target


#2 Features data (features)

#3 Distribution Features data (dist_features)

data.frame


Specifies the spatial distribution of conservation features across planning units. Each row corresponds to a combination of planning unit and feature.

  • pu
    integer id of a planning unit where the conservation feature listed on the same row occurs.

  • feature
    integer id of each conservation feature.

  • amount
    numeric amount of the feature in the planning unit. Set to 1 to work with presence/absence.

#pu feature amount
#1        3      1
#2        3      1
#3        3      1
#4        3      1
#5        3      1
#6        3      1

#3 Distribution Features data (dist_features)

#3 Distribution Features data (dist_features)

#4 Threats data (threats)

data.frame


Specifies the threats to consider in the optimization exercise. Each row corresponds to a different threats.

  • id
    integer unique identifier for each threat.

  • blm_actions (optional)
    numeric penalty of connectivity between threats. Default is 0.

  • name (optional)
    character `name for each conservation feature.

#>id    name blm_actions
#> 1 threat1           0
#> 2 threat2           0

#4 Threats data (threats)

#5 Distribution Threats data (dist_threats)

data.frame


specifies the spatial distribution of threats across planning units. Each row corresponds to a combination of planning unit and threat.

  • pu
    integer id of a planning unit where the threat listed on the same row occurs.

  • threat
    integer id of each conservation feature.

  • amount
    numeric amount of the threat in the planning unit. Set to 1 to work with presence/absence. Continuous amount values require that feature sensitivities to threats be established.

#>pu threat amount action_cost status
#> 8      2      1           2      0
#> 9      2      1           2      0
#>10      2      1           2      0
#>11      1      1           3      0
#>11      2      1           4      0
#>12      1      1           3      0

#5 Distribution Threats data (dist_threats)

data.frame


specifies the spatial distribution of threats across planning units. Each row corresponds to a combination of planning unit and threat.

  • action_cost
    numeric cost of an action to abate the threat in each planning unit.

  • status (optional)
    integer value that indicate if each planning unit should be available to be selected (0), locked-in (2) as part of the solution, or locked-out (3) and excluded from the solution. :::

#5 Distribution Threats data (dist_threats)

#6 Sensitivity data (sensitivity)

Optional

data.frame


specifies the sensitivity of each feature to each threat. Each row corresponds to a combination of feature and threat.

  • feature
    integer id of each conservation feature.

  • threat
    integer id of each threat.

  • a (optional)
    numeric the minimum intensity of the threat at which the features probability of persistence starts to decline.

  • b (optional)
    numeric the value of intensity of the threat over which the feature has a probability of persistence of 0.

#>feature threat
#>      1      1
#>      2      1
#>      3      1
#>      4      1
#>      1      2
#>      2      2

#6 Sensitivity data (sensitivity)

Optional

data.frame


specifies the sensitivity of each feature to each threat. Each row corresponds to a combination of feature and threat.

  • c (optional)
    numeric minimum probability of persistence of a features when a threat reaches its maximum intensity value.

  • d (optional)
    numeric maximum probability of persistence of a features in absence threats. :::

#6 Sensitivity data (sensitivity)

a, b, c and d

#6 Sensitivity data (sensitivity)

a, b, c and d

#7 Boundary data (boundary)

Optional

data.frame


Specifies the spatial relationship between pair of planning units. Each row corresponds to a combination of planning unit.

  • id1
    integer id of each planning unit.

  • id2
    integer id of each planning unit.

  • boundary
    numeric penalty applied in the objective function when only one of the planning units is present in the solution.

#>id1 id2 boundary
#>  1   1        0
#>  2   1        1
#>  3   1        2
#>  4   1        3
#>  5   1        4
#>  6   1        5

Workflow

Q3. Optimization objective?

Minimize Costs (minimizeCosts) vs Maximize Benefits (maximizeBenefits)


\(x_{ik}\) is the decisions variable that specifies whether an action to abate threat \(k\) in planning unit \(i\) has been selected \((1)\) or not \((0)\).

\(c_{ik}\) is the cost of the action to abate the threat \(k\) in the planning unit \(i\).

\(b_{is}(x)\) is the benefit of the feature \(s\) in the planning unit \(i\) (ranging between \(0\) and \(1\)).

\(t_s\) is the recovery target for feature \(s\).

\(f(x)\) is the function of connectivity penalty.

Q3. Optimization objective?

Minimize Costs (minimizeCosts) vs Maximize Benefits (maximizeBenefits)


\(x_{ik}\) is the decisions variable that specifies whether an action to abate threat \(k\) in planning unit \(i\) has been selected \((1)\) or not \((0)\).

\(c_{ik}\) is the cost of the action to abate the threat \(k\) in the planning unit \(i\).

\(b_{is}(x)\) is the benefit of the feature \(s\) in the planning unit \(i\) (ranging between \(0\) and \(1\)).

\(B\) is the budget available.

\(f(x)\) is the function of connectivity penalty.

Workflow

II) Create mathematical model

Create the model (problem())

Create an optimization model for the multi-action conservation planning problem, following the mathematical formulations used in Salgado-Rojas et al. (2020).


problem(
  x,
  model_type = "minimizeCosts",
  budget = 0,
  blm = 0
)

Workflow

III) Solve model

Solve model (solve())

Solves the optimization model associated with the multi-action conservation planning problem. This function is used to solve the mathematical model created by the problem() function.


solve(
  a,
  solver = "",
  gap_limit = 0,
  time_limit = .Machine$integer.max,
  solution_limit = FALSE,
  cores = 2,
  verbose = TRUE,
  name_output_file = "output",
  output_file = TRUE
)

Solve model (solve())

solver (gurobi, cplex, symphony)


Rsymphony

Solve model (solve())

Solve model (solve())

time limit (seconds)

solution_limit (First solution?)

verbose (Information on screen?)

Overview

Getting information of solutions

getActions()

Returns the spatial deployment of the actions for each planning unit of the corresponding solution.


getActions(x, format = "wide")

Output example,

##   solution_name pu 1 2 conservation connectivity
## 1           sol  1 0 0            0            0
## 2           sol  2 0 0            0            0
## 3           sol  3 0 0            0            0
## 4           sol  4 0 0            0            0
## 5           sol  5 0 0            0            0
## 6           sol  6 0 0            0            0

getSolutionBenefit()

Returns the total benefit induced by the corresponding solution. The total benefit is computed as the sum of the benefits obtained, for all features, across all the units in the planning area.


getSolutionBenefit(x, type = "total")

Output example,

##   solution_name feature benefit.conservation benefit.recovery benefit.total
## 1           sol       1                    0               11            11
## 2           sol       2                    0               16            16
## 3           sol       3                    0               10            10
## 4           sol       4                    0                9             9

getCost()

Provides the sum of costs to actions and monitoring applied in a solution.


getCost(x)

Output example,

##   solution_name monitoring threat_1 threat_2
## 1           sol         61       20       65

Outline

  • Features

  • Toy example

  • Overview

  • Actual case study: The Mitchell River, Australia

  • Shiny Application

Package installation and documentation

Package prioriactions can be found at CRAN, so it can be installed using:

install.packages("prioriactions")

Latest stable versions can be downloaded and installed from GitHub as follows (package remotes should be installed first):

if (!require(remotes)) install.packages("remotes")
remotes::install_github("prioriactions/prioriactions")

So, we load the prioriactions package.

# load package
library(prioriactions)

Toy Example

(Available in the Get Started vignette on the prioriactions website)


  • This example contains 100 planning units, 4 features and 2 threats.

  • The distribution of features and threats can be plotted on a grid of 10 x 10 units.

  • prioriactions contains this simulated example inside the setup files. You can extract it by the data() function.

#1 Planning units data

The monitoring cost values ranging from 1 to 10 and all status of 0 (not locked).

# load planning unit data from prioriactions
data(sim_pu_data) #To load simulated data

#head(sim_pu_data)

#1 Planning units data

A RasterLayer object can be used to present this spatial information. The pixel values correspond to the monitoring costs of each planning unit.

library(raster) #To plot rasters

r <- raster(ncol=10, nrow=10, xmn=0, xmx=10, ymn=0, ymx=10)
values(r) <- sim_pu_data$monitoring_cost
plot(r)

#2 Features data

Contains information about the features such as its id and targets (mandatory when minimizing costs).

# load features data from prioriactions
data(sim_features_data)

#head(sim_features_data)

#3 Features distribution data

Contains information on the spatial distribution of these features across planning units.

# load features data from prioriactions
data(sim_dist_features_data)

#head(sim_features_data)

#3 Features distribution data

To plot the spatial distribution of the first feature,

# load amount of features data
features <- reshape2::dcast(sim_dist_features_data, 
                            pu~feature,
                            value.var = "amount", 
                            fill = 0)

values(r) <- features$`1`
plot(r)

#3 Features distribution data

#4 Threats data

Provides information about the threats such as their id and name.

# load threats data from prioriactions
data(sim_threats_data)

#5 Threats distribution data

Provides information on the spatial distribution of these threats

# load threats data from prioriactions
data(sim_dist_threats_data)

#5 Threats distribution data

#6 Sensitivity data

Indicates which features is sensitive to what threat.

# load threats data from prioriactions
data(sim_sensitivity_data)

#7 Boundary data

Provides information on the spatial relationship between PU and they are presented in long format.

# load boundary data from prioriactions
data(sim_boundary_data)

Step 1: Initialize the problem

After having loaded our data, we will now create the data object through the inputData() function.

# create conservation problem
b <- inputData(pu = sim_pu_data, 
               features = sim_features_data, 
               dist_features = sim_dist_features_data, 
               threats = sim_threats_data, 
               dist_threats = sim_dist_threats_data, 
               sensitivity = sim_sensitivity_data, 
               boundary = sim_boundary_data)

# print problem
print(b)
## Data
##   planning units: data.frame (100 units)
##   monitoring costs:     min: 1, max: 10
##   features:       feature1, feature2, feature3, feature4 (4 features)
##   threats:        threat1, threat2 (2 threats)
##   action costs:   min: 1, max: 10

Step 1: Initialize the problem

It is advisable to determine the maximum benefit achievable for each conservation feature through the getPotentialBenefit() function

# get benefit information
getPotentialBenefit(b)

Step 2: Create the mathematical model

# create optimization problem
c <- problem(b, model_type = "minimizeCosts")
## Warning: The blm argument was set to 0, so the boundary data has no effect
## Warning: Some blm_actions argument were set to 0, so the boundary data has no
## effect for these cases
# print model
print(c)
## Optimization Problem
##   model sense: minimization
##   dimensions:  284, 396, 12.656 kB (nrow, ncol, size)
##   variables:   396

Step 3: Solve the model

To solve our model, we need an optimizer. In this case, we will use the Rsymphony library.

# solve optimization problem
d <- solve(c, solver = "symphony", verbose = TRUE, output_file = FALSE, cores = 2)

# print solution
print(d)
## Solution overview
##   name: sol
##   objective value: 738
##   gap:  0%
##   status:  Optimal solution (according to gap tolerance: 0)
##   runtime: 0.017 sec

Getting information

getActions()

# get actions from solution
actions <- getActions(d, format = "wide")

# print first six rows of data
# head(actions)

Getting information

getActions()

values(r) <- actions$`1`
plot(r)

More connectivity?

We go back to step 2 and repeat.


# create optimization problem
c2 <- problem(b, model_type = "minimizeCosts", blm = 10)
## Warning: Some blm_actions argument were set to 0, so the boundary data has no
## effect for these cases
# print model
print(c2)
## Optimization Problem
##   model sense: minimization
##   dimensions:  29984, 10296, 883.856 kB (nrow, ncol, size)
##   variables:   10296

More connectivity?

# solve optimization problem
d2 <- solve(c2, solver = "symphony", verbose = TRUE, output_file = FALSE, cores = 2)

# print solution
print(d2)
## Solution overview
##   name: sol
##   objective value: 863
##   gap:  0%
##   status:  Optimal solution (according to gap tolerance: 0)
##   runtime: 22.847 sec
# get actions from solution
actions2 <- getActions(d2, format = "wide")

values(r) <- actions2$`1`
plot(r)

Outline

  • Features

  • Overview

  • Toy example

  • Actual case study: The Mitchell River, Australia

  • Shiny Application

The Mitchell River, Australia

  • We divided the whole catchment (71,630 km2) into 2316 sites (i.e., sub-catchments)

  • We sourced the distribution of 45 fish species in the Mitchell river catchment as our conservation features.

  • We considered four major threats to freshwater fish species in the catchment: water buffalo (Bubalis bubalis), cane toad (Bufo marinus), river flow alteration (caused by infrastructure for water extractions and levee banks) and grazing land use.

The Mitchell River, Australia

Outline

  • Features

  • Overview

  • Toy example

  • Actual case study: The Mitchell River, Australia

  • Shiny Application

Shiny Application

Shiny Application

Thank you!

Questions?